%
% Author: Vasyl Mykhalchuk
%

function [ output_args ] = scalarField( anmFileName, frameNumber )
% Modify matlab path
addpath(genpath('../util'));
addpath(genpath('../io'));
addpath(genpath('../strain'));
addpath(genpath('../export'));
addpath(genpath('..'));

% Load animated mesh
anmFileNamePath = strcat('../../res/',anmFileName);
anmFileNamePath = strcat(anmFileNamePath, '.anm');
[geometryInfo, topologyInfo] = loadANM(anmFileNamePath);

%==========================================================================
dimensionsGeometry = size(geometryInfo);
dimensionsTopology = size(topologyInfo);

numberOfFaces = int32(dimensionsTopology(1));
numberOfFrames = int32(dimensionsGeometry(1));
numberOfVertices = int32(dimensionsGeometry(3) / 3);

if(frameNumber >= numberOfFrames)
    error('frameNumber should not exceed number of frames in animation');
end%if

%==========================================================================
% compute adjacency lists
adjList = createAdjacencyList(topologyInfo, numberOfVertices, numberOfFaces);
adjVertList = createVertexAdjacencyList(topologyInfo, numberOfVertices, numberOfFaces);

%==========================================================================

% time-varying scalar field (vertex based)
scalarFields = zeros(numberOfVertices, 1, 'double');
% time-varying scalar field (triangle based)
scalarFieldsT = zeros(numberOfFaces, 1, 'double');


for iFace = 1:numberOfFaces
    tr = topologyInfo(iFace,:);
    v1 = geometryInfo(1, 1, tr(1)*3 + 1: tr(1)*3 + 3);
    v2 = geometryInfo(1, 1, tr(2)*3 + 1: tr(2)*3 + 3);
    v3 = geometryInfo(1, 1, tr(3)*3 + 1: tr(3)*3 + 3);
    
    vv1 = geometryInfo(frameNumber, 1, tr(1)*3 + 1 : tr(1)*3 + 3);
    vv2 = geometryInfo(frameNumber, 1, tr(2)*3 + 1 : tr(2)*3 + 3);
    vv3 = geometryInfo(frameNumber, 1, tr(3)*3 + 1 : tr(3)*3 + 3);
    
    %flatten triangles 
    [v1f, v2f, v3f] = fitTrianglesToOnePlane(v1, v2, v3);
    [vv1f, vv2f, vv3f] = fitTrianglesToOnePlane(vv1, vv2, vv3);
    
    strains = tri2D_strain(v1f, v2f, v3f, vv1f, vv2f, vv3f);
    scalarFieldsT(iFace) = strains(1);
end%for

% compute vertex-based time-varying scalar field
deformationMax = 0.0;
for iVertex = 1:numberOfVertices
        i = int32(1);
        scalarValue = 0.0;
        adjTriangle = adjList(iVertex,i);
        while (adjTriangle ~= 0)
            tr = topologyInfo(adjTriangle,:);
            v1 = geometryInfo(frameNumber, 1, tr(1)*3 + 1 : tr(1)*3 + 3);
            v2 = geometryInfo(frameNumber, 1, tr(2)*3 + 1 : tr(2)*3 + 3);
            v3 = geometryInfo(frameNumber, 1, tr(3)*3 + 1 : tr(3)*3 + 3);
            trArea = 1; %triangleArea(v1,v2,v3);
            scalarValue = scalarValue + trArea * scalarFieldsT(adjTriangle);
            i = i + 1;
            adjTriangle = adjList(iVertex,i);
        end
        if scalarValue > deformationMax
            deformationMax = scalarValue;
        end%if
        scalarFields(iVertex) = scalarValue;
end%for
    
for i=1:numberOfVertices
        if deformationMax ~= 0
            scalarFields(i) = scalarFields(i) / deformationMax;
        end
end%for


targetScalarField = scalarFields;
[sizeX, sizeY] = size(targetScalarField);
%targetScalarField = reshape(targetScalarField, sizeY, sizeX);
%targetScalarField = avrgScalarField;
criticalPoints = extractCriticalPoints(targetScalarField, adjVertList);
%save('output/criticalPoints.txt', 'criticalPoints', '-ASCII');
colorList = createColorInfo(criticalPoints, targetScalarField);
vertexList = geometryInfo(1,1,1:numberOfVertices*3);
prepareMeshForExport(vertexList, topologyInfo,colorList,anmFileName);

end%function buildScalarFld





